# nice plots
library(lattice)

# open db and import data
library(RMySQL)
con = dbConnect(MySQL(), dbname="vmstats", user="openstack", pass="openstack")
vmstats = dbGetQuery(con, statement = 
                       'select vm,host,cpu,date_format(time, "%Y-%m-%d %H:%i:%S") as timestamp from vmstats where time < "2013-01-23" and time > "2013-01-18"')
vmstats$timestamp <- with (vmstats, as.POSIXct(strptime(timestamp, format="%Y-%m-%d %H:%M:%S")) )

# look into table
str(vmstats)
head(vmstats)

# get number of hosts
all.hosts <- unique(unlist(vmstats["host"]))
all.hosts

# get vm names
all.vms <- unique(unlist(vmstats['vm']))
all.vms

# write as csv
write.csv(vmstats,  file="vmstats.csv", row.names = FALSE)

# just for offline play
vmstats <- read.csv("vmstats.csv", stringsAsFactors = FALSE)
vmstats$timestamp <- with (vmstats, as.POSIXct(strptime(timestamp, format="%Y-%m-%d %H:%M:%S")) )

# get load per vm in a period
vmload <- vmstats[vmstats$vm == 'vm11',
                  c('timestamp', 'cpu')]
head(vmload)
plot(vmload$timestamp, vmload$cpu, type="p", cex=0.2, xlab="Time", ylab="CPU")


# smoothing
vmload$cpu <-filter(vmload$cpu, rep(1/60,60), side=1, circular=TRUE)
lines(vmload$timestamp, vmload$cpu, type='l', col='red')
abline(h=mean(vmload$cpu))
####

# all VMs in one plot
xyplot (cpu ~ timestamp | vm, groups=vm, data=vmstats, type=c("l", "g", "r"),
        layout=c(1, length(all.vms)), as.table=TRUE, xlab="Time", ylab="CPU")

# smoothed all VMs - not as nice as without 
xyplot (cpu ~ timestamp | vm, groups=vm, data=vmstats, type=c("r", "g", "spline"),
        layout=c(1, length(all.vms)), as.table=TRUE, xlab="Time", ylab="CPU")

####
# summary per host

host.load.sum <-function(host.name) {
  host.stats <- vmstats[vmstats$host == host.name, c('host', 'timestamp','cpu')]
  result <- aggregate(host.stats['cpu'], by=host.stats['timestamp'], FUN=sum, na.rm=TRUE)
  return(result)
}
s = host.load.sum('f18-3')

xyplot (cpu ~ timestamp, data=s, type=c('l'), as.table=TRUE, xlab="Time", ylab="CPU")
xyplot (cpu ~ timestamp, data=s, type=c("r", "g", "spline"), as.table=TRUE, xlab="Time", ylab="CPU")

##
# all hosts loads 
hosts.loads <- aggregate(vmstats['cpu'], by=c(vmstats['timestamp'], vmstats['host']), FUN=sum, na.rm=TRUE)
xyplot(cpu ~ timestamp | host , data=hosts.loads, type=c('g', 'l'), cex=0.1, layout=c(1,3), as.table=TRUE)
xyplot(cpu ~ timestamp | host , data=hosts.loads, type=c('g', 'r', 'spline'), cex=0.1, layout=c(1,3), as.table=TRUE)
###

# decomposition of time series
# f(t) = a + b * sin(t) + c * cos(t) with period = 24 hours
TIME <- unclass(vmload$timestamp)
COS <- cos(2 * pi * TIME / (24 * 60 * 60))
SIN <- sin(2 * pi * TIME / (24 * 60 * 60))
vmload.lm <- lm(vmload$cpu ~ COS + SIN )

summary(vmload.lm)
a <- coef(vmload.lm)[1]
b <- coef(vmload.lm)[2]
c <- coef(vmload.lm)[3]
my.f <- a + b * COS + c * SIN 

plot(vmload, type='p', cex=0.1)
lines(vmload$timestamp, my.f, col="red")
abline(h = mean(vmload$cpu), col="blue")

vm.amplitude = sqrt(b*b + c*c)
vm.amplitude
vm.phase = atan(b/c)
vm.phase

# check that function is really what we have found
my.f2 <- a + vm.amplitude * sin(2 * pi * unclass(vmload$timestamp) / (24 * 60 * 60) + vm.phase)
points(vmload$timestamp, my.f2, col="green", cex=0.1)

# corrected function if mean < amplitude
if (vm.amplitude > a) {
  a <- vm.amplitude
}
my.f2 <- a + vm.amplitude * sin(2 * pi * unclass(vmload$timestamp) / (24 * 60 * 60) + vm.phase)
points(vmload$timestamp, my.f2, col="red", cex=0.1)

# Re-arrangement is done one by one. We check what host has bigger peak loads 
# p = (sum(a) + sqrt( (sum(b))^2 +(sum(c))2 )
# Then we take all VMs of the host and try to find better home for them - calculating possible peaks

# vmloads

decompose.vmloads <- function (p.vmstats) {
  # get vm names
  vms <- unique(unlist(p.vmstats['vm']))
  a <- b <- c <- rep(0, length(vms))
  vmloads <- data.frame(a, b, c,stringsAsFactors = FALSE)
  rownames(vmloads) <- vms
  # get load per vm in a period
  for(vm in vms) {
    vm.stat <- vmstats[p.vmstats$vm == vm,
                       c('timestamp', 'cpu')]
    # smoothing
    vm.stat$cpu <-filter(vm.stat$cpu, rep(1/60,60), side=1, circular=TRUE)
    # decomposition of time series
    TIME <- unclass(vm.stat$timestamp)
    COS <- cos(2 * pi * TIME / (24 * 60 * 60))
    SIN <- sin(2 * pi * TIME / (24 * 60 * 60))
    vmstat.lm <- lm(vm.stat$cpu ~ COS + SIN )
    a <- coef(vmstat.lm)[1]
    b <- coef(vmstat.lm)[2]
    c <- coef(vmstat.lm)[3]
    ampl <- sqrt(b^2 + c^2)
    # move sin() up if amplitude > mean
    if(a < ampl) a = ampl;  
    vmloads[vm, 'a'] <- a
    vmloads[vm, 'b'] <- b
    vmloads[vm, 'c'] <- c
  }
  return (vmloads)
}

vms.coeffs <- decompose.vmloads(vmstats)
vms.coeffs


# check coefficients here, just to make nice pictures
TIME <- seq(min(vmstats$timestamp), max(vmstats$timestamp), by = 60*60)
timestamp <- rep(TIME, length(all.vms))
vm <- rep(all.vms, each=length(TIME))
cpu <- rep(0, length(vm))
time.int <- unclass(timestamp)

vm.stats <- data.frame(vm,timestamp, time.int, cpu, stringsAsFactors = FALSE)

cpu.calc <- function(x) {
  my.t = as.integer(x[3])
  my.vm = x[1]
  my.cpu = vms.coeffs[my.vm,'a'] + 
    vms.coeffs[my.vm,'b'] * cos(2 * pi * my.t / (24 * 60 * 60)) + 
    vms.coeffs[my.vm,'c'] * sin(2 * pi * my.t / (24 * 60 * 60))
  return(my.cpu)
}

vm.stats$cpu <- apply(vm.stats, 1, cpu.calc)

xyplot (cpu  ~ timestamp | vm, data=vm.stats, type=c("l", "g", "r"),
        as.table=TRUE, xlab="Time", ylab="CPU", layout=c(3,4))

# hosts peaks, calculated from vm coefficients
# p = (sum(a) + sqrt( (sum(b))^2 +(sum(c)) ^ 2) )
# host.vm is a table of hosts and corresponding vms

decompose.host.loads <- function (p.vms.coeffs, p.host.vms) {
  # get host names
  hosts <- unique(unlist(p.host.vms$host))
  a <- b <- c <- p <- rep(0, length(hosts))
  hosts.coeffs <- data.frame(a, b, c, p, stringsAsFactors = FALSE)
  rownames(hosts.coeffs) <- hosts
  for(host in hosts) {
    vms <- p.host.vms[p.host.vms$host == host,'vm']
    host.vms.coeffs <- p.vms.coeffs[vms,]
    hosts.coeffs[host, c('a', 'b', 'c')] <- c(sum(host.vms.coeffs$a),
                                              sum(host.vms.coeffs$b),
                                              sum(host.vms.coeffs$c))
    hosts.coeffs[host, 'p'] = hosts.coeffs[host, 'a'] + 
      sqrt(hosts.coeffs[host, 'b'] ^2 + hosts.coeffs[host, 'c'] ^2)
  }
  return (hosts.coeffs)
}

hosts.vms <- unique(vmstats[c('host', 'vm')])
hosts.coeffs <- decompose.host.loads ( vms.coeffs, hosts.vms)
hosts.coeffs

host.cpu.calc <- function(x) {
  my.t = as.integer(x[3])
  my.host = x[1]
  my.cpu = hosts.coeffs[my.host,'a'] + 
    hosts.coeffs[my.host,'b'] * cos(2 * pi * my.t / (24 * 60 * 60)) + 
    hosts.coeffs[my.host,'c'] * sin(2 * pi * my.t / (24 * 60 * 60))
  return(my.cpu)
}

hosts.to.optimize <- all.hosts

## repeat for all hosts ##
x <- function(){} # dummy

hosts.to.optimize
length(hosts.to.optimize)

# if in hosts to optimize only one is left - we are finished

# picture to check sin-like host loads
# TIME is calculated earlier
timestamp <- rep(TIME, length(hosts.to.optimize))
host <- rep(hosts.to.optimize, each=length(TIME))
cpu <- rep(0, length(hosts.to.optimize))
time.int <- unclass(timestamp)

hosts.stats <- data.frame(host,timestamp, time.int, cpu, stringsAsFactors = FALSE)
hosts.stats$cpu <- apply(hosts.stats, 1, host.cpu.calc)
hosts.stats

xyplot (cpu  ~ timestamp | host, data=hosts.stats, type=c("l", "g", "r"),
        as.table=TRUE, xlab="Time", ylab="CPU", layout=c(1,3))

### overloaded host

## plot current distribution of VMs
hosts.vms
vms.places<-table(hosts.vms, dnn=c("",""))
vms.places
levelplot(vms.places, col.regions=c("white", "blue"), border='black', colorkey=FALSE, xlab="", ylab="")

hosts.coeffs
hosts.coeffs.optimize <- hosts.coeffs[hosts.to.optimize,]
hosts.coeffs.optimize
overloaded.ind <- which (hosts.coeffs.optimize$p == 
                           max(hosts.coeffs.optimize$p),
                         arr.ind=TRUE)
overloaded.host <- rownames(hosts.coeffs.optimize)[overloaded.ind]
overloaded.host

# VMs to move
vms.to.move <- hosts.vms[hosts.vms$host == overloaded.host, 'vm']
vms.to.move
# minimize max peak
# start with the first vm
vm.to.move <- NA
vm.new.host <- NA
max.host.load <- max(hosts.coeffs.optimize$p)
for(vm in vms.to.move) {
  # if max load among hosts after vm move is less, replace vm.to.move
  # new load coeffs of host = old load coeffs + vm load coeffs
  other.hosts <- hosts.to.optimize[hosts.to.optimize != overloaded.host] 
  for(host in other.hosts) {
    new.load <- hosts.coeffs[host, 'a'] + vms.coeffs[vm, 'a'] +
      sqrt( (hosts.coeffs[host, 'b'] + vms.coeffs[vm, 'b']) ^ 2 +
              (hosts.coeffs[host, 'c'] + vms.coeffs[vm, 'c']) ^ 2 )
    if(new.load < max.host.load) {
      max.host.load <- new.load
      vm.to.move <- vm
      vm.new.host <- host
    }
  }
}

vm.to.move
vm.new.host
vm.old.host <- overloaded.host

# if we have no vm.to.move we cannot optimize further host with max load
if(is.na(vm.to.move )) {
  # remove host from optimiation attempts  
  hosts.to.optimize <- hosts.to.optimize[hosts.to.optimize != overloaded.host]
} 
hosts.to.optimize

# change hosts.coeffs, hosts.vms  
vm.a <- vms.coeffs[vm.to.move, 'a']
vm.b <- vms.coeffs[vm.to.move, 'b']
vm.c <- vms.coeffs[vm.to.move, 'c']

hosts.coeffs[vm.new.host, 'a'] <- hosts.coeffs[vm.new.host, 'a'] + vm.a
hosts.coeffs[vm.new.host, 'b'] <- hosts.coeffs[vm.new.host, 'b'] + vm.b
hosts.coeffs[vm.new.host, 'c'] <- hosts.coeffs[vm.new.host, 'c'] + vm.c
hosts.coeffs[vm.new.host, 'p'] <- hosts.coeffs[vm.new.host, 'a'] + 
  sqrt( hosts.coeffs[vm.new.host, 'b'] ^ 2 +
          hosts.coeffs[vm.new.host, 'c'] ^ 2 )

hosts.coeffs[vm.old.host, 'a'] <- hosts.coeffs[vm.old.host, 'a'] - vm.a
hosts.coeffs[vm.old.host, 'b'] <- hosts.coeffs[vm.old.host, 'b'] - vm.b
hosts.coeffs[vm.old.host, 'c'] <- hosts.coeffs[vm.old.host, 'c'] - vm.c
hosts.coeffs[vm.old.host, 'p'] <- hosts.coeffs[vm.old.host, 'a'] + 
  sqrt( hosts.coeffs[vm.old.host, 'b'] ^ 2 +
          hosts.coeffs[vm.old.host, 'c'] ^ 2)

hosts.vms[hosts.vms$vm == vm.to.move, 'host'] <- vm.new.host
hosts.vms
